Soft Sphere Packings at Finite Pressure but Unstable to Shear 
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When are athermal soft sphere packings jammed? Any experimentally relevant definition must 
at the very least require a jammed packing to resist shear. We demonstrate that widely used 
(numerical) protocols in which particles are compressed together, can and do produce packings 
which are unstable to shear — and that the probability of generating such packings reaches one 
near jamming. We introduce a new protocol that, by allowing the system to explore different box 
shapes as it equilibrates, generates truly jammed packings with strictly positive shear moduli G. For 
these packings, the scaling of the average of G is consistent with earlier results, while the probability 
distribution P(G) exhibits novel and rich scalings. 

PACS numbers: 05.70.Jk,05.10.-a,62.20.D- 
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Foams, emulsions, colloidal suspensions, granular me- 
dia and other particulate media undergo a jamming 
transition when their constituent particles are packed 
densely enough [IHZ]- This transition has been exten- 
sively studied in packings of deformable, athermal, fric- 
tionless spheres interacting through purely repulsive con- 
tact forces [BTI12). The limit where the particles just 
touch then plays the role of an unusual critical point, 
as a host of quantities, such as shear modulus, time and 
length scales, and contact number exhibit power law scal- 
ing with the distance to this critical point [8UT7]. 

Numerically created particle packings play a central 
role in many fields of physics, in particular jamming. In 
all numerical jamming studies we are aware of, packings 
are created by compressing a collection of particles, ei- 
ther by inflating the particles or shrinking the simula- 
tion box [8"Ml7|. It is then widely believed and tacitly 
assumed that, when compressed, the system simultane- 
ously develops a finite pressure, a finite yield threshold 
O [10] and a positive shear modulus G [8HT3] . Here we 
demonstrate that, to the contrary, algorithms that work 
solely by compression tend to produce packings that are 
unstable to shear, and hence have negative shear moduli. 
Nevertheless, such 'improperly jammed' packings possess 
a positive pressure P and a positive bulk modulus, and 
are in mechanical equilibrium — see Fig. la. 

In this Letter, we probe and explain this anomaly. The 
root problem is that compression only (CO) algorithms 
ignore the global shear degrees of freedom. We find that 
this results in a fraction of improperly jammed CO pack- 
ings which reaches one at the critical point. Hence, com- 
pression alone does not lead to jammed packings, and 
previous results on jamming have considered packings 
that, instead of being jammed, have been linearly unsta- 
ble to shear — in particular near jamming. 

Furthermore, we remedy this anomaly by introducing 
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FIG. 1: (Color online) (a) Example of a well-equilibrated CO 
packing of N = 32 particles which is unstable to shear (pres- 
sure P = 10 -2 , bulk modulus K ~ 0.385, contact number 
z w 4.26). (b) Illustration of the sinusoidal angular depen- 
dence of G on the principle direction of shear, 8, for three 
different packings at the same N and P — curve III corre- 
sponds to the packing shown in (a), and dashed lines indicate 
Gdc, the angular average of G. 



a shear stabilized (SS) packing algorithm that produces 
truly jammed packings with positive definite shear mod- 
uli [18] . and probe the probability distribution of G, un- 
covering novel scaling with distance to jamming and sys- 
tem size. 

Shear moduli in CO Packings — We have generated 2D 
packings of N soft harmonic bidisperse disks (with unit 
spring constant [TTJ) by a standard CO packing gener- 
ating algorithm, for pressures P ranging from 10~ 6 to 
10" 1 and 16 < N < 1024. Prior studies of the shear 
modulus have focused on ensemble averages at fixed dis- 
tance to the jamming point (P), typically for large N, 
and without reference to the angular dependence of G. 

As illustrated in Fig. lb, fluctuations and anisotropy 
are key: G varies sinusoidally with 9, and its angular- 
average, Gdc j varies substantially with realization. We 




FIG. 2: (Color online) Energy landscape where \r) denotes 
the particle degrees of freedom, and AX the box-shape. CO 
packings sit at a minimum of U with respect to |r); SS pack- 
ings sit at a minimum with respect to both \r) and AX. 



distinguish three types of packings. (I) Truly jammed 
packings for which G(6) > 0. (II) Improperly jammed 
packings for which G(9) < (III) Improperly jammed 
packings for which G{9) becomes negative over an in- 
terval in 9. We stress that all these packings are in a 
mechanical equilibrium and have a positive bulk modu- 
lus. 

It has been customary to measure G along a fixed di- 
rection [TUl [HI E0T - f2"4] . and the limited unstable range 
of type III packings, combined with the rare occurrence 
of type II packings, may explain why these instabilities 
have escaped attention to date. Since simulations often 
produce some "problematic" packings (for example due 
to issues with convergence) , packings of types II and III 
have likely been treated as "bad apples" and thrown out 
of the ensemble [2SJ [55] . 

Boundaries and Shear Stabilization - - Improperly 
jammed packings are not caused by numerical artifacts 
but stem from the essence of compression only (CO) al- 
gorithms. Consider the potential energy landscape as a 
function of the particle positions, \r), and shear defor- 
mations of the box, | AX) (Fig. 2). CO algorithms fix 
the unit cell and generate packings at a minimum of U 
with respect to \r). Residual shear stresses and shear 
moduli correspond to the first and second derivatives, 
respectively, of U along a strain direction AX — without 
permitting the strain degrees of freedom to equilibrate, 
both the residual stress and shear modulus are uncon- 
trolled. 

To create packings that are guaranteed to be stable 
against shear in all directions, we include shear deforma- 
tions of the box and search for local energy minima of 
U (Fig. 2) [27]. We combine standard conjugate gradi- 
ent techniques [26 with the FIRE algorithm [28 , which 
improves the speed by an order of magnitude, and also 
precisely control the pressure of the resulting packings. 
Since the energy is at a minimum with respect to the 
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FIG. 3: (Color online) Scatter plots of A m i n vs G for 50 pack- 
ings of N = 128 and P as indicated. Dots correspond to 
G(9 — 0), and blue (red) lines indicate the range of G(9) 
when the minimum of G(0) is positive (negative). The right 
bottom quadrant is empty: when A m i n > 0, G is positive def- 
inite, (a) SS packings, (b) CO packings at P = 10~ 2 . (c) 
CO packings at P — 10 -a — close to jamming, the fraction 
of improperly jammed CO packings grows dramatically. 



shear degrees of freedom, these packings have strictly 
positive values of G and exhibit zero residual shear stress 
[2"7] , unlike CO states. However, as a result of equilibrat- 
ing the strain degrees of freedom, the unit cell is no longer 
square. For example, starting from a CO packing (min- 
imum of U with respect to \r)), the box is deformed to 
find a minimum in the extended space spanned by \r) 
and the strain coordinates (Fig. 2). Such deformations 
are small for large systems |29j . 

A formal way of capturing the role of the boundaries is 
in terms of the stiffness matrices K° and K, where K° is 
the usual Hessian, while the "extended Hessian" K, in- 
troduced in a different context in Ref. [17], includes the 
dependence on the shear degrees of freedom — for details 
see the supplementary material. It can then be shown 
that G{9) is positive definite for all 8 if all eigenvalues of 
K are positive (excluding the trivial zero energy transla- 
tional modes). Defining A m in as the minimal eigenvalue 
of K, the sufficient condition for a packing to be stable 
against shear is A m j n > 0. In contrast, a positive spec- 
trum for the usual Hessian K° only guarantees stability 
in a box with fixed boundaries, but does not guarantee 
stability to all possible shear deformations (Fig. 1 and 2), 
contrary to the claim in Ref. [30] . 

Scatter plots of shear modulus and A m ; n for CO and SS 
ensembles shown in Fig. 3 confirm our picture: (i) All SS 
packings have positive A m j n and G. (ii) CO packings can 
have negative A m i n . Although many of these A m i„ < 
packings are stable when sheared along a fixed direction 
(dots correspond to 6 = 0), they almost always have 
negative G when sheared along other directions. 

Fraction of improperly jammed CO packings — What 
fraction of CO packings is unstable to shear? What gov- 
erns the scaling of this fraction? Fig. 4 shows that the 
probability that CO packings have shear directions along 
which G is negative, Pg<o , reaches one near jamming, 
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FIG. 4: (Color online) The fraction of CO packings unsta- 
ble to shear collapses when plotted as function of the excess 
number of contacts, NAz , where Az := ; 

z - 4 + 4/iV. 
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FIG. 5: (Color online) (a) Linear scaling of (G) with Az co for 
CO packings. The errorbars indicate the RMS fluctuations in 
G. (b) Linear scaling of (G) with Az ss +8/iV for SS packings 
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and that larger packings need lower pressures for these 
instabilities to become dominant. It is natural to expect 
that Pg<o would collapse when plotted as a function of 
L/l*, where I* is a characteristic length-scale which di- 
verges as 1/A.z near jamming, and where Az is the differ- 
ence between the contact number z and its value at the 
jamming point [HI IT21 IT51 I3lti3"3"] . Surprisingly, Fig. 4 
shows that the number of excess contacts ~ NAz, not 
the characteristic length scale /*, governs the fraction of 
improperly jammed packings — note that we have in- 
cluded a finite size correction to Az (see below). 

We conclude that the standard view of the jamming 
transition, in which rigidity is attained by simply com- 
pressing particles together [rTJHT^] . needs a revision: 
when the pressure is lowered in finite CO packings, more 
and more packings will become unstable to shear, leading 
to a blurring of the (un)jamming transition. We stress 
that many excess contacts are needed to avoid improp- 
erly jammed CO packings: for example, one needs of the 
order of a hundred excess contacts for Pg<o < 0.1. 

Scaling of Contact Number and G — Do the same scal- 
ing laws for, e.g., z or G [HJdl], govern both CO and SS 
packings? To answer this question, we have performed a 
finite size scaling analysis of both SS and CO packings: 
both the distance to jamming and the system size play a 
crucial role [34] . 

We first consider the contact number z [5HT21 155] . 
A packing is called isostatic when the number of con- 
straints, G, equals iVdof — No, the number of degrees of 
freedom A/^of minus the number of rigid body modes Nq. 
There is one constraint for each of the N c = Nz/2 force 
bearing contacts |36j . In two dimensions, No = 2, corre- 
sponding to two rigid body translations (rotation is in- 
compatible with periodic boundary conditions). Hence: 
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For CO states in two dimensions, iVdof = 2N (the particle 
displacements), so that z^° = 4 — A/N . For SS states 



the particle displacements are augmented by two shear 
degrees of freedom, leading to zf^ Q — 4. 

Is the isostatic bound reached at unjamming? We have 
found that both CO and SS packings have one contact in 
excess of their respective isostatic values when approach- 
ing the jamming point (see Suppl. Mat.). Goodrich et 
al. have argued that this extra contact reflects the re- 
quirement that jammed states have positive bulk modu- 
lus, which puts an additional constraint on the box size 

We now turn our attention to the scaling of G, and 
first investigate the scaling of the angle-averaged shear 
modulus (G) in ensembles of finite sized CO and SS pack- 
ings. In Fig. 5a we show that in the CO ensemble, (G) 
is proportional to z — z^®, consistent with prior results 
QUI EH ED EH 133- In Fig. 5b we show that in the SS 
ensemble, the average shear modulus is proportional to 
z — (z?q — 8/N). So, although the SS shear modulus is 
also linear in z, its vanishing point extrapolates to a state 
with four contacts less than the isostatic state. We note 
that in both ensembles (G) is of order 1/N in the zero 
pressure limit. 

The amount of scatter in (G) observed in our new CO 
packings is surprisingly large. We note that previous 
work did not consider the value of G over all angles and 
discarded negative values of G, which leads to a smaller 
scatter [351 US]- Recent work by Goodrichef al. shows 
that this scatter can be further suppressed by using ex- 
ceptionally accurate equilibration and larger ensembles 
[37] ■ Nevertheless, the observation that SS data exhibits 
far lower scatter than CO data, while both packings were 
obtained with the same numerical accuracy, suggests that 
remnants of the unstable modes present in the CO en- 
semble hinder accurate equilibration. 

With few exceptions [TUl fT5| fT6 | [38ti42] . studies of jam- 
ming have focused on ensemble averages. Here we con- 
sider the probability distribution P(G) for both ensem- 
bles, sampling both and realizations. Fig. 6a illustrates 
that for CO packings, P(G) often peaks at negative G, 
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FIG. 6: (Color online) (a) The probability distributions for G 
of CO and SS packings differ qualitatively, (b) Scaling of the 
variance ((G — (G)) 2 } for SS packings reveals novel scaling, 
(c) P(G/{G)) shows a systematic variation with LAS. 



and can possess an extended tail towards negative G. In 
contrast, for SS packings, G is strictly positive, and the 
peak of P{G) is always at finite G. 

For SS packings, the distributions P(G) are well- 
behaved; however, there is no single parameter scaling. 
For brevity of notation, we define z = 4 — 8/N, so that 
(G) ~ z — S = AS. Our data shows that the variance Oq 
scales roughly linear with Az/L (Fig. 6b). The scalings 
of the average and variance of G suggest that distribu- 
tions of P(G/(G)) that have equal values of LAS might 
collapse. Fig. 6c shows that grouping P(G/(G}) by LAS 
captures the main trends: for large LAS, the distribu- 
tion P(G/(G)) is clearly peaked away from zero, but for 
lower values of LAS becomes more skewed and wider. 
We note that the scaling of P(G < 0) for CO packings 
suggest that finite size scaling corrections for P(G) differ 
between CO and SS packings, and an important question 
for the future is to probe these differences [43] . 

Discussion — Improperly jammed CO packings domi- 
nate in the critical, near jamming regime, whereas pack- 
ings made by a shear stabilized algorithm are strictly 
jammed: boundary conditions play a crucial role in con- 
trolling the rigidity of packings, in particular close to 
jamming. In most experimental procedures, the creation 
history is richer than homogeneously inflating particles, 
and involves the motion of boundaries and shear [TJ |3j- 
[H] - - how does this relate to our scenario? First, 
we note that in contrast to the 'shear jammed packing' 
of Bi et al. [H], our CO and SS packings only exhibit 
small contact anisotropies that vanish as 1/%/iV [Mj, and 
that CO packings exhibit similarly weak anisotropies in 
the contact forces — we use shear to stabilize, rather 
than jam. Second, we note that the strong anisotropy 
of G that we observe is reminiscent of fragility as intro- 
duced by Cates et al., although usually fragile states are 
defined as having no resistance to shear in certain direc- 
tions (i.e., G = 0), while here we have G<0. Moreover, 
such fragility typically arises due to the shear history of 
the system [HI [35]. Nevertheless, it is conceivable that 
protocols that do not explicitly perform shear stabiliza- 



tion initially yield improperly jammed states, which then 
relax until they reach a fragile state. 
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an apparent force felt by particles involved in boundary- 
crossing contacts when the lattice vectors are distorted. 
It comprises the first 2N components of K\T), where 
|T) = |{0}2jv,7, 0). The quadratic term in the change in 
potential energy, which governs linear stability, is then 
MJ/V = {l/2V){q\K\q) = (1/2)G(<9) 7 2 , where V is the 
volume, so that 



G(9) = (q\K\q)h 2 V 



(3) 



From Eq. (pi) it immediately follows that G{9) is pos- 
itive definite for all 9 if all eigenvalues of K are positive 
(excluding the trivial zero energy translational modes). 

We finally note that Eq. ^ requires the extended Hes- 
sian K, which includes shear degrees of freedom [T7]. To 
relate G to the usual Hessian K°, we follow Ref. [21] and 
write G7 2 /2V = W a — W na , where W a > is the work 
done in affinely displacing the particles and the box. The 
actual particle displacements \u) have a non-affine con- 
tribution |u na ) = \u) — \u a ) that reduces the deformation 
energy by W na = (u na \K a \u na )/2. W na > if the spec- 
trum of K° is non-negative, but G can still be negative if 
W na > Ma- — hence there is no simple relation between 
the sign of G and the eigenvalues of the usual Hessian 
K°. 

Contact Numbers for P — > — Based on our count- 
ing, the isostatic number of contacts, C{ so , equals 2A^ — 2 
for the CO ensemble, and precisely 2N for the SS ensem- 
ble. In Fig. 7 we show our results for the excess contact 
numbers C + , defined as the difference between the actual 
number of contacts C and the respective isostatic value 
Ci S0 . For both the CO and SS ensembles, the number 
of excess contacts reaches one in the limit of vanishing 
pressure: so for 2D CO ensembles, the number of contacts 
reaches 2N — 1, and for 2D SS ensembles, the number of 
contacts reaches 2A^ + 1. 



Supp Material 

Extended Hessian — A packing's linear response to 
shear can be expressed as a matrix equation. Let us 



introduce \u) 
and 



IK,«niLi>. \i) 



iK,ur}£ 1>7 , 



d 2 U 



Kmn du m du n 



and K mn = 



d 2 U 



dq n ' 



(2) 



evaluated at the coordinates corresponding to a packing. 
K° is the usual Hessian or stiffness matrix, while the 
"extended Hessian" K, introduced in Ref. [IT], includes 
the dependence on 7 and 9. The response to imposed 
strain is the solution to K°\u) = \Fp), where \Fp) is 
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FIG. 7: (Color online) Scaling of the excess contact number 
of contacts, C+ = C - C iso for (a) the CO and (b) the SS 
ensemble. 



